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Abstract. Velocity profiles of ground state lines of HbD" 1 ", HC 18 + and N2H + , observed previously with the CSO 
and IRAM 30m telescopes, are modeled with a Monte Carlo radiative transfer program to study the temperature, 
density and velocity structure of the pre-stellar core LDN 1544. The IFD -1 " line is double-peaked like that of the 
other ions, but previous models that fit the HC ls O + and N2H + profiles are found not to fit the H2D 4 " data. 
Matching the H2D + observations requires at least three modifications to the model at small radii: (1) the density 
profile must continue to rise inward and not flatten off toward the center; (2) the gas temperature must be nearly 
constant and not drop inwards significantly; (3) the infall velocity must increase inward, in a fashion intermediate 
between 'quasi-static' (ambipolar diffusion) and 'fully dynamic' (Larson- Penston) collapse. The C ls O emission 
indicates a chemical age of <0.1 Myr. The effects of a flattened structure and rotation on the line profiles are 
shown to be unimportant, at least on the scales probed by single-dish telescopes. Alternatively, the H2D + profile 
is affected by absorption in the outer layers of the core, if gas motions in these layers are sufficiently small. 
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1. Introduction 

The formation of low-mass stars occurs in dense cores 
inside molecular clouds, although the actual mechanism 
by which these cores accumulate mate rial is still dis- 
puted (ambipolar diffu sion vs. turbulence: IShu et all 2000; 
IHartmann et al.ll200ljh In any case, through the loss of 
turbulent and/or magnetic support, dense cores contract 
towards a 'critical state', after which gravitational collapse 
starts and infall occurs onto a central object, a 'protostar'. 
Observers distinguish such protostars in Class objects, 
where the circumstellar envelope is much more massive 
than the star and accretion rates are high, and Class I 
objects, where star and envelope have similar masses and 
accretion rates have dropped (see lAndre et alJl200ol for a 
review). During these stages, most («99%) of the angular 
momentum of the core is carried away by bipolar out- 
flows, the formation of binary stars, and magnetic fields. 
The remaining 1% leads to the formation of the circum- 
stellar disks commonly o bserved in subseque nt stages of 
pre-mam sequence stars l|Mundv et all l200C|h This sce- 
nario is convincing, but many fundamental questions re- 
main unanswered regarding the stages before formation of 
a central luminous object, the so-called pre-stellar cores. 
Does the collapse start from a state of dynamical equi- 



librium? When does a rotating structure form, and how 
important are the effects of flattening? 

Detailed kinematical studi es have been carrie d out for 
the Class I object LDN 1489 faogerheiidell200ll) and the 
Class object IRAM 04191 l|Belloche et alJl2002|) . These 
authors used the HCO + and CS molecules, respectively, 
to trace the velocity field. However, at the low tempera- 
tures ( 5 10 K) of pre-stellar cores, most neutral molecules 
freeze out onto dust grains. In particular, the CO and 
CS molecules freeze out at n(H2) Z 10 5 cm -3 , so that 
these species (as well as HCO + ) trace the outer part s 
of the cores l(Bacmann et alJ 120031 iTafalla et ahll2004h . 
The more volatile N2 molecule remains abundant at these 
densities, making N2H + and N2D + bet ter tracers of pre- 
stellar core nuclei than C O and CS fe.g.. lBergin fc LangeJ 
ll997tlCrapsi et al.ll2004h . However, at the very centers of 
the cores, where Z 10 6 cm -3 , even N 2 freezes out, caus- 
ing N?H + and N?D + to disappear from th e gas phase 
llBergin et al J 120021: iBelloche fc Andrei l2004h . One of the 
few molecules that do not suffer from this depletion effect 
is H2D+. This paper explores the use of H2D+ as kine- 
matic tracer of the centers of pre-stellar cores. We focus on 
the 'prototypical' pre-stellar core LDN 1544 in the Taurus 
molecular cloud (d=140 pc) which has been well charac- 
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Table 1. Centroids and widths of the two velocity compo- 
nents. Numbers in brackets denote uncertainties in units 
of the last decimal. 



Line 


^LSR 


AV bs 






km s _1 


km s _1 


km s _1 


H 2 D+ (lio-ln) 


7.07(3) 


0.29(7) 


0.28-0.34 




7.38(3) 


0.27(4) 




HC 18 + (1-0) 


7.04(1) 


0.18(3) 


0.10-0.12 




7.28(1) 


0.23(3) 




N 2 H+ (1-0, FiF=10-ll) 


7.08(1) 


0.19(1) 


0.11-0.13 




7.33(1) 


0.20(2) 




a Thermal line width at T k 


n =7 and 10 K. 





terized by observation s and models. Preliminary results of 
this work appeared in lVan der Tak et al.l l)2004|) . 

The structure of this paper is as follows. After the 
data and the radiative program have been introduced 
(§§Ell3i we test models of quasi-static collapse mediated 
by ambipolar diffusion, which fit previous observations of 
LDN 1544, in spherical (§ 0} and axial (§ EJl symmetry. 
More general spherically symmetric free-fall models are 
discussed in §0 Sect.0tests the adopted molecular abun- 
dance structure of the core against chemical models. The 
possibility of dynamic rather than quasi-static collapse is 
considered in § The sensitivity of the results to the 
adopted collisional cross sections for H2D+ is examined in 
§EH Sect. EHhsts the conclusions and future prospects of 
this work. 



2. Data 

The line profile of H2D + observe d toward the 'dust peak' 
of LDN 1544 in October 2002 l|Caselli et alJ 120031) sug- 
gested a double-peaked shape but was of poor quality with 
low signal to noise. Therefore the source was included in 
our December 2003 CSO observations, which will be pre- 
sented by Caselli et al. (in prep.). The left panel of Figure^ 
shows the combined data at the central position, which 
has an rms noise level of T m b=0.11 K per 0.04 km s _1 
channel. In this spectrum, the H2D + line profile is clearly 
double-peaked, like th e profiles of the HC 18 Q+ and N 2 H+ 
1-0 lines, observed bv lCaselli et al.l 1 2002al) and shown in 
the central and right panels of the figure. To quantify the 
similarity between these lines, we fitted two Gaussians to 
the profiles and found that each line is well described by 
two thermal components at Tki n =7-10 K, separated by 
w0.26 km s" 1 (TableCJ. 

The H2D + line has also been detected at positions 
20" offset from the dust peak, where the line is only 
half as strong, suggesting a centrally peaked H2D+ abun- 
dance l)Caselli et all2003l) . Figure shows th e line profiles 
of H2D 4 - and of HC 18 0+ and N 2 H+ (from I Caselli et all 
2002a). as averages of observations taken at 20" North, 
South, East and West offsets. For H 2 D + and N 2 H + , nei- 
ther these 'average offset' spectra appear double-peaked, 
nor the individual spectra. In contrast, the average 



HC ls O + spectrum at the offset position is clearly double- 
peaked. 

The H 2 D+ line profile towards the dust peak could be 
affected by absorption in the outer parts o f the core, as 
seen in the J =1 -0 lines of CS and HC O+ l|Tafalla et alJ 
Il998h and N 2 H+ liWilliams et al.lll999l). Absorption is not 
seen for DCO+ 2-1 l)Caselli et al.ll2002ajl but this may be 
an excitation effect. Analogously, the low H2D + intensity 
and lack of strong central absorption at the offset position 
may be an excitation effect. To clarify these points, de- 
tailed radiative transfer models have been run for H2D+, 
as described in the following sections. 



3. Model setup 

The line profiles of H 2 D+, HC ls O+ and N 2 H+ toward 
LDN 154 4 have been modeled with the Mon te Carlo pro- 
gram bv iHogerheiide fc van der Takl (j^OOQ) 1 , which can 
treat spherical and axisymmetric geometries. For HC 18 + 
and N?H + , we use m olecular data from the database by 
ISchoier et all l|2005h 2 . For H 2 D+, collisional rate coeffi- 
cients have not been calc ulated, so we use scaled radiative 
rates IjBlack et alJll990l) . which probably have an uncer- 
tainty of a factor of 3-10. Section investigates the sensi- 
tivity of the model results to the adopted rate coefficients. 

Modeling of the N2H + J=1^0 line concentrates on 
the lowest frequency FiF=10^11 hypcrfinc component, 
which has an optical depth of only 3.7% of the total tran- 
sition. We treat this line as if it were an isotopic version of 
N2H+: molecular excitation and line formation are calcu- 
lated with the N2H+ abundance reduced by 27. Similarly, 
to model the HC 18 + J=l— >0 line, an oxygen isotope 
ratio of 16 O/ 18 O=500 was adopted. We do not model 
H 13 CO + spectra beca use the line profile is comp licated 
by hyperfine structure (ISchmid-Burgk et alJ l2004). 



4. Spherically symmetric models 

The spherical models have inner and outer radii of 
80 and 13000 AU, on a 100-point grid, logarithmically 
spaced. Figure E3 shows the adopted temperatu re and den- 
sity structure, taken from iGalli et al.l l(2002l) . as well as 
the adopted molecular abundance profiles. For HCO + 
and N9 H + , we use abundance profiles from ICaselli et ail 
(2002b) (Model 3 with a low ionization rate and a 
high sticking coefficient, which fits their data best), and 
assumed zero abundance inside r=2500 AU. The dis- 
tribution of HC 18 + has a central hole because its 
chemical precursor, CO, freezes out onto dust grains at 
n(H2)>10 5 cm -3 . The N2H+ distribution continues fur- 
ther inward and has a hole at smaller radii, where N 2 
freezes out. The assum ption of an NoH + 'ho le' of 2500 AU 
was already needed bv ICaselli et alJ l)2003|) to reproduce 
the large H2D+ abundance at the center of LDN 1544. 
For H2D+, we initially assumed a central abundance of 



1 http : //www.mpif r-bonn.mpg . de/ staf f /f vandertak/ratran/ 

2 http : //www. strw.leidenuniv.nl/~moldata/ 
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Fig. 1. Line profiles of H 2 D + (left), HC 18 + (middle) and N 2 H + (right; only the Fx F=l, 0-1,1 hyperfine component) 
observed towards the ' dust peak' of LDN 154 4. Superposed are synthetic profiles for the best-fit static model (dotted), 
and infall models after ICiolek fe Basul l|2000|) for times 2.660 (dashed) and 2.684 Myr (dash-dotted). The line width 
(static: &£>=0.15 km s _1 ; infall: 6d=0.05 km s _1 ) is set to match the H2D + data. 



1x10 9 , dr opping by a factor of 5 at a radius of 20" 
(2800 AU) llVan der Tak et a Jl2004|) as suggested by the 



observations of ICaijellT'eT'aTT 

For the velocity field, we first adopt static models, 
with zero velocity at all radii, and second models with 
'infalP velocity fields. There are several competing the- 
ories for the velocity field of collapse onto a pr otostar. 
Here we consider models bv lCiolek fc Basul l|2000|) . which 
assume that protostellar collapse is quasi-static and reg- 
ulated by magnetic pressure. Specifically, we model their 
simulations 't3' and 't5' which correspond to times of 2.660 
and 2.684 Myr after the start of core collapse. At earlier 
times, the infall speed is very low, while later times are 
implausible due to the short time spent in these phases. 
We consider these specific infall models because they are 
known to reprodu ce observations of H CO + and N2H + to- 
ward LDN 1544 (jCaselli et all 120024 - Sections El and El 
treat more general collapse models. Note that the ve- 
locity field in the simulations by Ciolek & Basu is not 
spherically but axially symmetric with an axis ratio of 
«0.3 (see their Figure 4). Adopting these velocity fields 
in spherical symmetry is therefore a simplification; how- 
ever, the assumed profiles are quite similar to those of 



collapsing: B onnor-Ebert spheres, as recently shown by 
iMverj (|2005^ . Section will present fully consistent, two- 
dimensional models with the velocity field of Ciolek & 
Basu. For the turbulent broadening, Doppler parameters 
bo between 0.05 and 0.25 km s _1 were tried. Smaller val- 
ues of bo are overwhelmed by thermal broadening; larger 
values do not fit the data. Variations of bjj with radius are 
considered in Sectional 



Figure ^ shows that the total widths of the H2D" 1 " , 
HC ls O + and N2H+ ground-state lines can be matched 
using either velocity field. The best-fit static model 
has &£>=0.15 km s _1 , while &c=0.05 km s _1 gives the 
best fits with the infall models. However, the observed 
double-peaked line shape rules out the static model, 
both at the dust peak and at the 20" offset position. 
Using &£>=0.05 km s _1 the infall models can reproduce 
the HC 18 + and N2H+ line profiles at both positions 
(Figure 0). The data do not allow to decide between the 
't3' and 't5' solutions. However, none of these models re- 
produces the double-peaked H2D + profile seen at the cen- 
tral position. 
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Fig. 2. As Figure [2 at 20 arcseconds offset (average NSEW) 



One reason for the lack of a central dip in the syn- 
thetic H^D -1- spectrum may be that the model has H2D" 1 " 
abundant at th e core center, where infall velocities are 
low according to lCiolek fc BasvJ l|2000h . Therefore we have 
tested the possibility that H2D+ has a shell- type distri- 
bution like HC 18 0+ and N 2 H+. The H 2 molecule is too 
light to freeze out on dust grains and cause depletion of 
H2D+, but instead, conversion into D2H + and D^~ in the 
gas phase (e.g., H2D + +HD^D2H + +H2) may cause a cen- 
tral Ii2D + hole. Formation of multiply deuterated H^j" is 
expected at t he temperature and density of the cente r 
of LDN 1544 ^Roberts et alJl200fl IWalmslev et alJl2004h . 
Support for these theories comes from the recent detec- 
tion of tow ards the pre-stellar core LDN 1689N 
(jVastel el 

We have run models with the H2D + abundance set 
to zero inside inner radii of 500-2500 AU (in steps of 
500 AU) , but none of these models gives a double-peaked 
H2JD+ line profile. The problem is that the lCiolek fc Basul 
(2000) models predict infall velocities of <0.1 km s _1 
at the radii where H^D" 1 " is abundant, much less than 
the observed peak separation of w0.26 km s _1 (Table QJ. 
Increasing the outer radius of the Ii2D + shell from 
3000 AU to 6000 AU gives too strong line intensities at 
the offset position, but does produce double-peaked H^D" 1 " 
profiles at the center. However, the two peaks are not 



^equally bright as observed, but the blueshifted peak is 
sstwice as strong as the redshifted one, a phenomenon 
known as 'infall asymmetry'. Therefore it seems that the 
H2D + line shape cannot be explained with the physical 
and geometrical structure adopted so far. Before chang- 
ing the adopted physical structure of the core, we explore 
deviations from the hitherto assumed spherical symmetry. 



5. Axisymmetric models 

Maps of dust continuum and molecular line emis- 
sion of LDN 1544 show a flattened s tructur e 
l|Ward-Thomnson et alJ Il999t ICaselli et all l2002ah . 
For example, the 1300 /jm dust continuum, which should 
be representative of the surface density distribution, has 
a major axis of «2'.5 (21000 AU), an axis ratio of rs1.7, 
and a position angle of 45° for the major axis. Therefore 
we explore the effects of flattened geometry on the line 
profiles. 

We have constructed two-dimensional radiative trans- 
fer models with inner and outer radii of 70 and 20000 AU, 
on a 30x30 grid, logarithmically spaced in both R and 
z. The density structure is as described bv lCiolek fc BasiJ 
(200CJ): an 'infall disk' with a flaring angle of 15° and an 
inclination of 74°. We adopt a time of 2.660 Myr, which 
gives the best fit to millimeter-wave dust continuum ob- 
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Fig. 3. Radial profiles of dust and gas temperature, and the densities of H 2 , H 2 D + (x2xl0 9 ), HCO + (xlO 12 ) and 
N 2 H+ (xlO 10 ) in our model of LDN 1544. 



servations l|Ciolek fc Basul2000h . The radial infall velocity 
is zero at cloud center and edge, and peaks at 0.14 km s _1 
at a radius of 0.03 pc. For the gas and dust tempera- 
tures, we tried a uniform value of 10 K (as in t he model of 
Ciolek and Basu) and the radial distribution bv lGalli et alJ 
(2002]) (assumed vertically isothermal). The abundance 
profil es of N?H+, HC 18 Q+ and DCO+ were again taken 
from ICaselli et all l|2002b(l . Model 3, and for H 2 D+, we 
used the profile described in § 

These models reproduce the observed sp atial distribu- 
tion o f the DCO+ 2-1 and N 2 H+ 1-0 lines l|Caselli et all 
l2002a|) . and also predict the correct line shape for HC 18 + 
and N 2 H + . However, they predict a single-peaked H 2 D + 
line, which is not observed (Fig.HJ. 

Another possible cause of the double peaked H 2 D + 
line profile is rotation al motion of the cent ral part of 
the core. For example, iBelloche et al found the 

Class object IRAM 04191 to rotate at an angular rate 
of £1 ^3.9xl0~ 13 s _1 inside a 'centrifugal radius' rc of 
3500 AU. We have adde d rotation to the velocity field 
bv ICiolek & Basil l)2000l) . while keeping the density and 
temperature structure from the previous models. Rotation 
rates between fi=lxl0 -13 and2xl0 -12 s _1 were tried, in- 
side rc=3000 and 6000 AU. Larger values of rc are ruled 
out by the lack of a velocity shift of the H 2 D + profiles at 
the offset positions, and also by the inter f erome tric obser- 
vations of N 2 H+ 1-0 bv IWilliams et alJ l|l999t) . We find 
that the H 2 D + line profile is hardly changed from the 
non-rotating case for r c =3000 AU. For r c =6000 AU, the 
H 2 D + line is considerably broadened, and the total width 
of the observed H 2 D + profile limits f2 to <5xl0~ 13 s _1 . 
None of these models predict double-peaked line shapes. 
In general, both flattening and rotation appear unimpor- 




-1 -0.5 0.5 1 

Velocity (km s ! ) 



Fig. 4. Results of 2D models with T kin =10 K (dotted line) 
and with the temperature distribution from Galli et al 
(dashed line) , superposed on the observations (histogram) . 

tant for LDN 1544 on the scales probed by single-dish 
observations. 

6. Alternative temperature and density structures 

The reason that the models in §§0] and |S] do 

not indicate H 2 D + self- absorption may be that the 
H 2 D + excitation temperature is almost constant in the 
range of radii ( < 3000 AU) where H 2 D + is abundant 
(H 2 D + /H 2 = 10~ 9 ). In particular, the models from § 0] 
with Galli's temperature distribution have T ex =5.8-5.6 K 
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inside r=3000 AU (Fig. [5] top) , while the isothermal mod- 
els of §E]have T cx =8.8-8.2 K at these radii. 



Galli et al model 




10 

Radius (AU) 
Evans & Young model 




10 10 10 

Radius (AU) 
Isothermal model (T = 8 K) 




10 10 
Radius (AU) 

Fig. 5. Excitation temperature profiles of H2D + lio~lii 
for the physical structures by Galli et al (top), Young 
et al (middle), and for the isothermal power-law model 
(bottom). Within the radius of high H2D + abundance 
(~3000 AU), T ex is practically constant in the Galli case, 
but has a gradient in the other cases. 



In search of an excitation temperature gradient, 
we went back to spherically symmetric models and 
adopted the density and dust temperature distributions 



bv lEvans et al 



bv lYoung et al 



2001), with gas temperatures calculated 



(2004). In these models for LDN 1544, the 



density continues to rise toward the center, rather than 
flatten ing off as in the model of IWard-Thompson et all 
l)l99fll) . As a consequence, T cx drops from 6.9 to 5.5 K at 
radii of 500-2500 AU (Fig. El middle). 

Combi ning this physical stru cture with velocity field 
't3' from ICiolek & Basi] l|200(t . the H 2 D+ profile re- 
mains single peaked. This result holds also when adopt- 
ing constant infall speeds up to 0.2 km s _1 . Larger in- 
fall speeds are inconsistent with the HC 18 + and N2H + 
data. Even an inward increasing velocity profile, v — 
Vo(r/ro) ~ ' 5 with ro=10 4 AU and vq—0.1 km s _1 , fails 
to produce a double peak (Fig. EJ. Adopting ro=10 3 AU 
and fn=0.2 k m s , which mimics the 't5' velocity field 
dCasellil EorB does not change this result appreciably. 
Furthermore, applying central holes in the H2D + abun- 
dance distribution, with radii of 500-1500 AU, changes 
the total intensity of the line, but not its shape. 

To prevent T ex from dropping off at the very center, 
we took the more drastic step of forcing a uniform Tkin, 
while keeping the density power law structure by Evans 
et al. Values of T kin =6 K, 8 K, 10 K and 12 K were tried. 
Fig. (bottom) shows the T ox profile for T kin =8 K. The 
H2D + profiles for these models show a classic 'infall asym- 
metry' for abundances > 10 -9 , or a broad peak for lower 
abundances. The observed profile with two peaks of equal 
strength is not reproduced. Leaving the dust temperature 
unchanged or setting it equal to T^ n makes no difference 
to these results. 




-0.5 







0.5 



Velocity (km s ) 

Fig. 6. Results of models with a power-law density dis- 
tribution and Tkin=8 K (dotted line) and with the gas 
temperature calculated by Young et al (dashed line), su- 
perposed on the observations (histogram). 
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7. Self-consistent chemical model 

The rate of formation of H 2 D + by the reaction of 
with HD depends mainly on the cosmic-ray ionization rate 
which is not expected to change within the pre-stellar core. 
However, the main destruction channel of H2D+ is the re- 
action with CO, the abundance of which is expected to 
decrease significantly at low temperatures and high den- 
sities due to freeze-out on dust grains. The models so far 
take depletion into account in a simplified way by consid- 
ering 'jumpy' H2D + abundance profiles. This section con- 
siders models that compute the abundances of and its 
isotopomers self-consistently across the LDN 1544 core. 
For that we adapted the model of Hj~ chemistry in young 
protop lanetary disks developed bv ICeccarelli fc Dominikl 
( 2005) . These models ignore the effect of dynamics on the 
chemistry. Briefly, the model computes the abundances of 
H^", H2D+, D2H" 1 " and D^" by solving the chemical network 
for these species, which is mostly governed by the CO de- 
pletion and the cosmic-ray i onization rate l|Roberts et alJ 
120031 IWalmslev et all l2004t see the detailed discussion 
in lCeccarelli fc Dominikll2005h . The model uses the 'fast' 
rate for the H3 + HD reaction ijRoberts et alJEoTK^ and 
920 K for the binding e nergy of CO on the grain surface 
IjBergin fe Langerll997|) . The papers by Roberts et al. and 
Walmsley et al. discuss the uncertainties of the other rel- 
evant chemical reation rates such as the dissociative re- 
combination of H2D+, D2H + and D^. The role of cosmic 
rays is twofold: on the one hand, they set the ionization 
degree in the gas, and therefore the total abundance of H^ 
and its isotopomers in the region where CO is depleted. 
The larger the ionization degree, the larger the absolute 
abundance of H^ in all its isotopic forms. Second, the 
cosmic rays regulate the desorption of CO molecules from 
icy grain mantles back into the gas phase. The model in- 
cludes 'spot heating' of grains by cosmic ray impacts, but 
no non-therm al desorption me chanisms which are much 
less effective l(Shen et al.ll2004|) . The larger the cosmic- 
ray flux, the higher the CO gas-phase abundance, and the 
lower the abundance of the deuterated forms of H^~. 

The model assumes that at time t—0, all CO is in 
the gas phase, so that the CO abundance is equal to 
the canonical value 9.5 x 10~ 5 with respect to H 2 (e.g., 
iFrerking et al.lll982T) . As time passes, the CO molecules 
condense out onto the grain mantles and disappear from 
the gas phase. The time scale of the CO depletion is set 
by the density and temperature, which are taken equal 
to those computed in the models of Galli et al. and 
Evans et al., previously described. The result is a zone 
where the increased CO depletion leads to an enhanced 
H2D + /H j" ratio, namely an incr eased H 2 D + abundanc e 
(see also ICaselli et all I2OO3I and ICeccarelli et alJl2004|) . 
Our computed H2D + , D2H + and D t abundances coin- 
cide with those calcul ated by iRoberts et alJ ll200.il) and 
I Walmslev et alJ (|2004fl when the same densities are con- 
sidered. 

Figure H shows an example of the computed H2D" 1 ", 
D2H + and D^ abundance profile across LDN 1544. Also 
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Fig. 7. Calculated abundances of H 2 D + , D 2 H + and D^" 
for an age of 0.1 Myr, a cosmic-ray ioniz ation rate of 
3xl0~ 17 s _1 and the physical structure by lEvans et alJ 
l|200lh . The dotted lines indicate the radii equivalent to 
one and two CSO beams, for comparison with observa- 
tions. 



shown is the electron fraction, which equals the total 
ionization. In these highly shielded regions, most charge 
comes from molecules, not from metals, which are als o 
expected to be highly depleted llCaselli et alJ l2002bl) . 
Note that in the very inner region, the H 2 D + /H^" and 



D2H + /Hg" abundance ratios have a 'plateau', whereas 
the D^/H^" abundance ratio keeps increasing. The re- 
sult is a 'hole' in the H 2 D + and D 2 H + abundances, 
because Djj" becomes the charge carrier when the CO 
depletion is larger than ~30. The result that T>% is 
the dominant i on at the cent e rs of pre-s tellar cores was 
found before bv lRoberts et al.1 l)2003(l and lWalmslev et all 
ll2004h; adopting the tem perature and density structure 
from lYoung et alJ l|2004h for LDN 1544, it happens at 
i?^ 3000 AU. The H2D + line profiles computed from these 
models depend on two major parameters: the cosmic- 
ray ionization rate and the time (for the CO deple- 
tion). We ran models varying both parameters for the 
two physical structures of the core, the Galli et al. 
and Evans et al. structure respectively. We considered 
ages of 0.05, 0.1 and 1.0 Myr, and fixed the cosmic-ray 
ionization rate at the canonical v alue of 3xl0 -17 s _1 
(IVan der Tak fc van Dishoeckfe oOOl. The model does no t 
consider any shielding of cosmic rays, as in lCasellil l|2003f) . 

Since the central CO-free zone grows in time, the 
strengths of isotopic CO li nes constrain the chem ical age. 
With the velocity field of ICiolek fc Basul |)2000(K we cal- 
culate an intensity for the C 18 1-0 line in the 22" beam 
of the IRAM 30m telescope of 4.3 K for t=0.05 Myr, 
3.5 K for 0.1 Myr and 1.5 K for 1 . Myr . The ob- 
served strength of «5 K l|Caselli et al.lEo02a|) thus sug- 
gests an age of <0.1 Myr for LDN 1544. It is not clear 
if these ages are compatible with our assumption of a 
static clo ud. Models including both chemistry and dynam- 
ics (e.g., iRodgers fc Charnlevl EooS lAikawa et al ~l l2005h 
are needed to address this issue. 
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To model the H2D 4 " line profile, we assume an or- 
tho/para ratio of unity as appropriate for low tempera- 
tures and high densities, although the exact amount of 
molecular depletion may make a factor of ~2 difference 
l|Pagani et alJl992HWalmslev et all2004|) . Figure El shows 
the H2D + line profiles calculated for a constant infall ve- 
locity of 0.2 km s _1 and for the Ciolek & Basu velocity 
field. The synthetic profiles are double peaked, but show 
strong infall asymmetry while two peaks of equal strength 
are observed. The strong self-absorption makes the emis- 
sion «50% weaker than observed. At the offset position, 
the emission is also predicted to be asymmetric, and also 
«50% weaker than observed. 



E-i 



H 2 D+ llO _1 ll 






-0.5 







0.5 
-Is 



Velocity (km s ) 

Fig. 8. Results of models with a power-law density distri- 
bution and the gas temperature calculated by Young et 
al, the effect of depletion included for a chemical age of 
0.1 Myr for either a constant infall velocity (dotted line) or 
the Ciolek & Basu velocity field (dashed line), superposed 
on the observations (histogram). 



These models seem to overestimate the depletion 
of CO and the abundance of H2D+. Chemical ages 
<0.05 Myr do not seem plausible, but the age is uncer- 
tain through the dust opacity coefficient, which makes the 
density uncertain by a factor of two. Alternatively the gas 
temperature (which controls the emission) is higher than 
the dust temperature (which controls the depletion). The 
binding energy of CO to the grain surface may also have 
been overestimated, so that CO does not stick to the sur- 
face as well as assumed. We have not considered the re- 
duced cosmic-ray ionization rate which previ ous studies of 
the LDN 1544 core indica te in this source ijCaselli et alJ 
l2002bl iTafalla et alJEool . since the CO depletion would 
be even higher, implying a worse match to the data. 



8. Alternative kinematics 

The infall model bv ICiolek & Basul l(2000F ) is a magneti- 
cally mediated, 'quasi-static' collapse with relatively slow 
infall velocities. The collapse may instead be 'dynamic' 
which means that the gas never reaches equilibrium. The 
most extrem e case is the Larson-Penston collapse (e.g., 
IZhoulll992h where the velocity reaches many times the 
sound speed. We modeled the H 2 D + line with such a ve- 
locity profile, but find that emission is predicted at much 
higher velocities than observed. 

Hydrodynamic simulations bv lHunte'rl l)l977(l show the 
existence of a continuum of solutions between the extreme 
static and dynamic cases. The solutions for the density p 
at time t and radius r are 



AnGpr 2 



~2at 



and the velocity u (directed inward) follows 



it 
a 



where a is the sound speed, and mo is a dimensionless 
parameter. The Larson-Penston collapse has mo = 46.9 
while mo = 0.975 for the quasi-static (Shu) collapse. 

As examples of infall models 'intermediate' between 
quasi-static and fully dynamic collapse, we have run mod- 
els for Hunter's cases lib and lid which have mo = 2.577 
and m = 1.138 respectively. The gas and dust temper- 
atures in these models are equal at 10 K, corresponding 
to a sound speed of 0.19 km s _1 . Fig. [§] presents syn- 
thetic line profiles for Model lid. Results are shown for 
the case of a constant H2D + abundance and of an inner 
hole in the H2D 4 " distribution. These models have peaks 
at the observed velocities, but with 'infall asymmetry': the 
blueshifted peak is stronger than the redshifted one. The 
results of Model lib are similar, but with this asymmetry 
more pronounced. 

9. Effect of the collisional rate coefficient 

Another potentially important parameter for the H 2 D + 
line formation is the collisional rate coefficient (CRC). To 
quantify the sensitivity of our model results on the CRC, 
we have run models with the CRC increased and decreased 
by factors of 3 and 10. These models use the power-law 
density pro file fromlEvans et alJ l(200ll) . the temperature 
profile from Young et al.l (120041) and the H2D + abundance 
profile from lCeccarelli fc Dominikl l|2005|) . Since the previ- 
ous sections indicate that the details of the infall velocity 
field play a minor role for our models, we simplify the 
kinematics to a 'step function' in bp, the turbulent line 
width. The observations suggest that bo is larger in the 
(emitting) central part of the core than in the (absorbing) 
outer layers, so we adopt 6d=0.1 km s _1 at radii <10 3 AU 
and bo—0 outside this radius. The model does not include 
any systematic radial motions, but the increased 6d at the 
center is presumed to arise in infall. This setup implies the 
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Fig. 9. Results of models with a power-law density dis- 
tribution, Tki n =10 K, Hunter's velocity profile 'lid', and 
with the H2D + abundance constant (dotted line) and with 
an inner hole (dashed line), superposed on the observa- 
tions (histogram). 



existence of infall motions, while the results do not depend 
on the specific infall model. Note that thermal broadening 
is added to the turbulence at all radii. 

Figure shows the results of models with the CRC 
increased by facto rs of 3 and 10 above the value from 
Blac k et al I lll990h . The line strength is seen to scale with 
the CRC; this trend continues for the models with the 
CRC decreased below 'standard' which are not shown. 
The shape of the H2D + profile is seen not to depend on the 
CRC: both the emission from the core and the absorption 
in the outer layers scale with the CRC. For the adopted 
model parameters, the best fit to the observed H2D+ pro- 
file would be obtained by increasing the CRC by a factor 
of ~5. However, this number depends critically on several 
parameters, particularly n(H2), 7ki n and the H2D + abun- 
dance. This is why we refrain from recommending a value 
for the CRC at this point. Even though n(H.2) and T^ n can 
be obtained from other observations (as we have done), the 
H2D + abundance cannot. Therefore, for the time being, 
H2D + abundances can only be determined to factors of 
~10 accuracy. We recommend theoretical calculations of 
the CRC of H2D+ to improve this situation. 

The fact that the simplified kinematics give a better fit 
to the observed H2D" 1 " profile than the models in Sectional 
implies that the infall motions toward LDN 1544 are 
very small, in good a greement with the conclusion from 
Wil liams et alJ (|1999l) that the infall speed is much smaller 
than the thermal, rotational and gravitational speeds in 
this source. The observed intensity at the 20" offset posi- 
tion of T m b=0.4 K is well reproduced by the same model 
that fits the central spectrum. The model predicts a «20% 



absorption at the line center, which the present data do 
not rule out. However, better spectra at offset positions 
are needed to confirm this prediction. 
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Fig. 10. Results of models with the temperature den- 
sity distributions from Young et al, a 'jump' in bo at 
i?=10 3 AU, and with the H 2 D + collisional rate coefficient 
at the standard value and increased by factors of 3 and 10 
(curves from bottom to top), superposed on the observa- 
tions (histogram). 



10. Conclusions 

Radiative transfer modeling of line profiles of H2D" 1 ", 
HC 18 0+ and N 2 H+ in the pre-stellar core LDN 1544 
shows that previous descriptions of the temperature, den- 
sity and velocity structure, that reproduce HC ls O + and 
N2H+ data, do not fit H2D + data, which probe smaller 
radii. At least two modifications are required at these 
radii, which affect H 2 D+, but not HC 18 0+ or N 2 H+. 

First, to make the excitation temperature of H 2 D + in- 
cr ease inward, the density cann ot flatten off as proposed 
bv lWard-Thompson et al] (1999), but must continue to in- 
crease as in the models bv lEvans et alJ ((2001) . In addition, 
only very small inward decreases of the gas temperature 
are allowed. Observations of NH3 toward a sample of pre- 
stellar cores also suggest a constant kinetic temperature 
with radius (Tafalla et al. I2002L EoO^) . The same result is 
obtained from ob servations of CO iso topes toward the pre- 
stellar core B 68 llBergin et al.ll2005]) . suggesting that cur- 
rent models of pre-stellar cores either underestimate the 
dust temperature because of changes in the dust opacity, 
or overestimate the gas-grain thermal coupling because 
of grain growth. However, the NH3 and CO observations 
probe larger radii than H2D + , and interferometric obser- 
vations of NH3 would be valuable. 
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Second, the infall velocity needs to increase with ra- 
dius down to small radii. If cloud collapse is quasi-static 
and mediated by ambipolar diffusion, the central region 
supported by thermal pressure, where the infall speed 
drops to zero, must be smaller than in the models by 
ICiolek fe Basul ((2000) . Alternatively, the collapse may be 
'mildly' d ynamic, indicated by values of 1-3 for th e param- 
eter m l|Hunterlll977t iFoster fe Chevalierlll993]l . Highly 
supersonic infall velocities, as well as large-scale rotation, 
are ruled out. 

Alternatively, the central dip on the H 2 D + profile is 
due to absorption in the outer parts of the core, as seen 
for CS, HCO + and N2H + . However, such absorption is 
not seen for DCO + 2-1, nor in H2D + at offset positions. 
Also, to produce an absorption as narrow as observed, the 
outer layers of LDN 1544 must have essentially zero infall 
and turbulent moti ons, consistent with e vidence from CS, 
HCO+ and N 2 H+ llWilliams et alJll999lh To test this 'ab- 
sorption' hypothesis, sensitive mapping of the H2D 4 " line 
and of DCO + 1-0 are needed. 

In the future, making full use of H2D + as kinematic 
probe requires maps at higher sensitivity and spatial and 
spectral resolution than the CSO can offer. Other single- 
dish telescopes (JCMT, APEX) can fulfill some of these 
requirements, but not all. Therefore, for a real break- 
through, interferometric observations will be crucial. The 
SMA will enable case studies, but only ALMA will have 
the sensitivity to probe many pre-stellar cores within rea- 
sonable times. The interpretation of the data will require 
good observational constraints of the gas temperature and 
density. In addition, theoretical calculations of the colli- 
sional rate coefficients of H2D+ are urgently needed. 

If the current observations cannot be matched with any 
of these modifications, then we speculate that the doubled- 
peaked H2D 4 " line profile in LDN 1544 is caused by the 
presence of two protostellar condensations, orbiting each 
other, in the core center. We may be witnessing the birth 
of a binary protostellar system. Once again, interferometer 
observations are required to test this hypothesis. 
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